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ABSTRACT 

We study the growth of density perturbations in a universe with unstable dark 
matter particles. The mass {niu) range 30 eV < rriu < lOkeV with lifetimes {td) in 
the range 10'' sec < < 10^^ sec are considered. We calculate the COBE normalized 
matter power spectrum for these models. We find that it is possible to construct 
models consistent with observations for masses rriu > 50 eV by adjusting so as to 
keep the quantity m^(keV)td(yi') constant at a value around 100. For m,y < IkeV 
the power spectrum has extra power at small scales which could result in an early 
epoch of galaxy formation. We do not find any value of td which gives a viable model 
in the mass range rrii, < 50 eV. We also consider the implications of radiatively 
decaying neutrinos — models in which a small fraction B <^ 1 of neutrinos decay into 
photons-which could possibly ionize the intergalactic medium (IGM) at high redshift. 
We show that the parameter space of decaying particles which satisfies the IGM 
observations does not give viable models of structure formation. 



Subject headings: large-scale structure of the universe — elementary 
particles:clustering — dark matter:intergalactic medium 
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1. 



Introduction. 



To understand the formation of large scale structures in the universe is one of the most 
challenging problem in standard big bang cosmology. In recent years an overwhelming amount 
of evidence has accumulated which suggests that the formation of structures can be understood 
within the framework of the gravitational instability picture in which the structures form through 
the growth of small, scale-invariant, adiabatic perturbations created in the very early universe 
by some process like inflation (see Liddle & Lyth 1993 and references therein). The consistency 
of the CMBR anisotropies measured by the COBE-DMR (Smoot et al. 1992; Bennett et al. 
1996) — which probe the matter fluctuations at recombination — with the present distribution of 
matter on large scales supports this view and suggests that the present universe is dominated 
by cold dark matter (CDM) particles. However the standard CDM model with = I, h = .5, 
and n = 1 normalized to COBE observations predicts neither the correct amplitude nor the 
correct shape of the power spectrum of density fluctuation at small scales (Efstathiou, Sutherland 
& Maddox 1990; Efstathiou, Bond, & White 1992; Peacock & Dodds 1994). The r.m.s. mass 
fluctuation in randomly placed spheres of radius 8h~^ Mpc (denoted by ag) provides a sensitive 
probe of the power spectrum at scales around k = 0.2/iMpc~^, and various studies show that the 
observed abundance of rich clusters of galaxies at the present epoch is consistent with cig ~ 0.5-0.8 
(Henry & Arnaud 1991; White et al. 1993; Viana & Liddle 1996; Bond & Myers 1996 Eke et al. 
1996; Pen 1996; Borgani et al. 1997; Carlberg et al. 1997), whereas the standard CDM model 
normalized to the four-year COBE data predicts erg = 1.22 (Bunn &; White 1996). In addition, 
the standard CDM model is also inconsistent with the shape of the power spectrum inferred from 
various galaxy surveys (Baugh & Efstathiou 1994; Lin at. al. 1996). 

Within the framework of CDM-like models, the power spectrum of density fluctuation can be 
analyzed using a single parameter (Bardeen at. al. 1986; Bond 1996) 



Here ilro = Pro/PcO, where p^o is the present energy density in all the relativistic species and 
Pco is the present critical density, and Jl^o is the contribution from the relativistic species to the 
present density parameter. Similarly, Q^o is the contribution to the present density parameter 
from the matter which can 'effectively' clump at small scales. While the standard CDM model 
predicts F = 0.5, observations require that 0.22 < L < 0.29 (Peacock & Dodds 1994), and several 
variants of the standard CDM model have been proposed to overcome this discrepancy. There 
are essentially two ways to decrease F: 1. by decreasing the amount of matter which can clump 
at small scales (i.e., by decreasing ^Imo) or 2. by increasing the radiation content of the universe. 
Models like ACDM (Kofman et al. 1993; Liddle et al. 1996a; Stompor et al. 1995 ), HCDM 
(Klypin et al. 1993) and oCDM (Liddle et al. 1996b; Yamamoto & Bunn 1996) all involve a 
decrease in the value of ^mO- Models with decaying dark matter particle achieve the required 
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decrease in T by increasing Qj-o (Bardeen, Bond, & Efstathiou 1987; Bond & Efstathiou 1991; 
White, Gelmini, & Silk 1995; McNally & Peacock 1996). In this paper we study the decaying 
neutrino model in some detail. 

Decaying neutrinos models are characterized by two free parameters — rrii,, the mass of the 
decaying neutrino and t^, the lifetime of the particle. In the early stages of the evolution when the 
temperature of the universe is much higher than rrii,, the energy density of the massive neutrinos 
is the same as the contribution from a massless species. As the universe expands the massive 
neutrinos become non-relativistic and the energy density in this component starts increasing 
relative to the contribution from the massless neutrinos. This gives rise to an era when the density 
of the universe is dominated by massive neutrinos, after which these particles decay and pump 
their rest energy into relativistic particles. Once the neutrinos decay the evolution of the universe 
is similar to the standard CDM model except that the energy density in relativistic particles is 
higher. This extra energy in the relativistic decay products gives the requisite increase in QrO (Eq. 
||) which acts to lower the value of T, and by suitably choosing rrii, and td it is possible to get 
results which are in agreement with the observed large scale structure. 

Decaying neutrino models have been investigated by earlier authors (Bond &: Efstathiou 1991; 
White, Gelmini, &: Silk 1995) who find that at the length-scales which are relevant for large scale 
structure formation {k <^ O.lMpc^^) the power spectrum is essentially of a CDM-like form with 
an effective T parameter which has a dependence on mu and td- This dependence captures the 
effect of the enhancement of i^ro^ a-nd in the mass range {mu > IkeV) which they have considered 
this is the only process which affects the shape of the power spectrum in the relevant range of 
k. However, if one considers neutrinos with smaller masses and larger lifetimes, there are other 
physical processes which become important and the power spectrum starts to exhibit features 
which cannot be described by just changing the T for a CDM-like power spectrum. For large 
masses the viable models are restricted to have short lifetimes for otherwise the energy density 
of the decay product is too high and the resultant value of F turns out to be too low. In such 
models the massive neutrino dominated era ends much before any of the relevant length-scales 
enter the horizon. For low neutrino masses the viable models require larger lifetimes, and for 
many such models the size of the horizon in the massive neutrino dominated era is comparable to 
the length-scales which are relevant for large scale structure formation. The main effect on the 
power spectrum is that the modes which enter the horizon in the massive neutrino dominated era 
have more power compared to a CDM-like power spectrum and a CDM-like fit based on adjusting 
F does not work for these models. 

In this paper we study in some detail the power spectrum for a large class of decaying 
neutrino models in the mass range 30 eV < rrii^ < lOkeV. In addition to the large masses 
considered in earlier works (Bond & Efstathiou 1991; White, Gelmini, & Silk 1995) we have also 
considered models with low masses and large lifetimes, the upper limit on the lifetime coming 
from the restriction that the massive neutrino has to decay before the present epoch. We have not 
considered models where the massive neutrinos survive until the present epoch. While for large 
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masses the free-streaming of the massive neutrinos has no effect on large scale structure formation, 
this process has to be taken into account when dealing with neutrinos with low masses and large 
lifetimes. Sub-horizon scale perturbations in relativistic collision-lcss particles are wiped out due 
to the free-streaming of the particles. The massive neutrinos behave like relativistic particles as 
long as their momentum is much larger than the rest mass, and sub-horizon scale perturbations 
in the massive neutrino component start growing only after the particles become non-relativistic. 
The earlier works dealt with the perturbations in the massive neutrino component using the 
hydrodynamic equations which cannot capture the effects of free-streaming, and hence they were 
restricted to large masses. In this paper we have used the collision-lcss Boltzmann equation which 
allows us to follow the free-streaming of the massive neutrinos and hence we have been able to 
study models with low masses and long lifetimes. For a detailed discussion of free-streaming the 
readers is referred to Bond &; Szalay (1983). The relevant equations for following the growth of 
perturbations in a decaying massive neutrino scenario are presented in the Appendix. 

The power spectrum for the various decaying neutrino models have been normalized so that 
the theoretically predicted CMBR anisotropics are consistent with the COBE-DMR observations. 
For some of the decaying neutrino models the normalization differs significantly from the standard 
CDM model and this happens because of a large contribution from the 'integrated Sachs- Wolfe 
effect'. For a discussion of the 'integrated Sachs-Wolfe effect' (Sachs k. Wolfe 1967) in the 
context of some other large scale structure formation models the reader is referred to Kofman &; 
Starobinsky (1985) and Sugiyama &; Gouda (1992). 

The motivation for studying cosmological models with decaying massive neutrinos has been 
twofold. As discussed earlier, decaying neutrinos provide possible scenarios for large scale structure 
formation. Decaying neutrinos are also interesting in the context of another very important 
problem in cosmology — to understand the ionization of the intergalactic medium (IGM) at high 
redshifts (Sciama 1990; Sethi 1997). The ionization state of the IGM at high redshifts is inferred 
from Gunn-Peterson tests for neutral hydrogen, neutral helium, and singly ionized helium, and 
also the proximity effect (see eg. Miralda-Escude k. Ostriker 1990; Giroux and Shapiro 1996; 
Sethi & Nath 1996), and it is not clear if the conventional sources of photoionization can serve the 
purpose of ionizing the IGM. Radiatively decaying neutrino models where a small fraction S <C 1 
of neutrinos decay into photons have been proposed as a possible means of photoionizing the 
IGM (Sethi 1997). In this paper we have investigated if any of the radiatively decaying neutrino 
models which can explain the ionization state of the IGM can also predict a large scale structure 
formation scenario which is compatible with observations. 

For the radiatively decaying neutrinos, in addition to the two parameters and td, the 
branching ratio for decay into photons B is another free parameter. However, since B this 
affects only the ionization state of the IGM. The decay photons are dynamically unimportant 
and they do not affect the evolution of the dark matter perturbations. Thus, for the purpose of 
calculating the matter power spectrum we need consider only mu and t^, as the relevant parameters 
and the value of B is of no consequence. 
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Prom the point of view of particle physics models, we require: (1) td <C to, to being the 
present age of the universe and (2) for radiatively decaying neutrinos, radiative lifetime (= t^/B) 
such that the intergalactic medium can be ionized. The latter was studied in Sethi (1997) and 
such models are possible within the framework of left-right symmetric models or models involving 
a charged Higgs scalar (see Pukugita and Yanagida 1995 for a recent review). To satisfy the 
first condition, we need to construct models where all the neutrinos have decayed by the present 
epoch. One of the models in which it can be implemented is a non-minimal Majoran model. In 
this model the dominant decay mode of the neutrino is a massless neutrino and a Majoran, and it 
is easy to get lO"^ < td < 10^^ sec for lOkeV > > 100 eV (Gelmini, Nussinov & Peccei 1992). 
The models we study are quite distinct from another class of radiatively decaying neutrino models 
suggested to solve several galactic and extragalactic observations by Sciama (1990). The dominant 
mode of decay in Sciama's neutrinos is radiative with radiative lifetimes of ~ 10^^ sec with masses 
nil, — 30 eV. As these neutrinos are essentially stable over the age of the universe, structure 
formation with these neutrinos closely resembles HDM models (for more details see Sciama 1994) 
and they are not considered here. 

Finally, we briefly outline the organization of the paper. In Section 2. we qualitatively discuss 
the physical effects that decide the power spectrum in a decaying neutrino model. The detailed 
equations and the computational method we use to numerically calculate the power spectrum are 
presented in the Appendix. Section 3. contains the results of our computations. We analyze the 
general features of the power spectrum for different decaying neutrino models and viable models 
are chosen on the basis of comparison to observations. In section 4. we present the summary 
and a discussion of the results. Appendix A contains the equations which govern the evolution of 
perturbations in the metric and in the various constituents of a decaying neutrino model. We work 
in the synchronous gauge and in our analysis we treat only the scalar perturbations. Appendix B 
contains a discussion of the initial conditions that we have used for the perturbations. Appendix C 
contains the equations which we have used to calculate the CMBR anisotropy and normalize the 
power spectrum. In Appendix D we briefly describe the scheme that we have used to numerically 
compute the power spectrum. 



2. Physical effects in the decaying neutrino model. 

In this section we qualitatively discuss some of the physical processes which shape the power 
spectrum in a decaying neutrino model. The detailed equations are presented in the Appendix. 

The first ingredient in any such model is the power spectrum of initial perturbations on scales 
larger than the horizon. Here we assume that the initial perturbations have been produced by 
some viable model of inflation and we use a scale invariant Harrison-Zcl'dovich power spectrum 
of the form Pi{k) = Ak. As the universe evolves the horizon expands, and the shape of the power 
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spectrum on sub-horizon scales gets modified by various astrophysical processes. The evolution of 

the power spectrum in the linear regime is expressed using the transfer function T{k) which relates 
the final power spectrum Pf{k) to the initial power spectrum i.e. Pf{k) = T{k)Pi{k). Below we 
discuss the various physical processes that decide the shape of the matter transfer function in a 
universe with decaying neutrinos. 

In all the models that we have studied CDM particles are the most dominant component of 
the present universe and it is the evolution of perturbations in this component which is of primary 
interest. The CDM perturbations are coupled to the other constituents of the universe through 
the gravitational potential which at any epoch is largely due to the most dominant component of 
the universe at that epoch. So we first discuss which of the various components dominates the 
universe in the different stages of its expansion and later discuss how this decides the shape of the 
transfer function. 

In a spatially flat universe the evolution of the scale factor is governed by the equation 



^a{v) = y ^P{v)a\v) = Ho^) . (2) 



where ri(t) = fQ{l/a{t')) dt' is the conformal time and i?o(= 100/ikm/s/Mpc) is the present value 
of the Hubble parameter, p{r)) is the total energy density of the universe (for decaying neutrino 
models it is contributed by the photons, 2 masslcss neutrinos, CDM particles, massive neutrinos, 
and neutrino decay product) and uj{t]) = p{ri)a^(7])/ pco, where pco is the present value of the 
critical density. The variable ^(r/) has been chosen so that at present it coincides with the value 
of the density parameter Qq = 1. 

We first consider a universe with only two kinds of constituents 1. relativistic particles 
(referred to as radiation) which at present contribute flrO to the density parameter, and 2. 
pressureless massive particles (referred to as matter) which at present contribute fi^o to the 
density parameter. The radiation make a constant contribution LOr{ri) = QrO while the contribution 

from matter uj^niv) = (^{v)^mO increases as the universe expands, and the universe proceeds from a 
radiation dominated era to a matter dominated era as it expands. The transition between the two 
regimes occurs when the matter and radiation make equal contributions to a'(??), and the value of 
the conformal time at the epoch of matter-radiation equality is given by 



r/eq = 2(V2-l)-f^. (3) 

In the radiation dominated era the Jeans length is of the order of the horizon (~ ci]), and 
perturbations in all components grow (i.e. 5{k,r]) oc r/^, where the mode k = 27r/A) when they 
are on scales larger than the horizon (i.e. they satisfy kcrj <C 1). The growth stops if the modes 
enter the horizon (i.e. k ~ 7r/(cr])) in the radiation dominated era. Once the universe gets 
matter dominated, matter perturbations on all scales grow in the same way {S{k,ri) oc rj"^). Using 
^eq = 7'"/(c??eq) to denote the mode which enters the horizon at the epoch of matter-radiation 
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equality, we can say that all modes with k < keq entered the horizon in the matter dominated era 
and they all grow by the same factor during the course of the evolution. The modes with k > keq 
enter the horizon in the radiation dominated era and as a consequence they experience a period 
of stagnation when they don't grow and the amplitude of these perturbations is suppressed by 
the factor (keq/k)"^ relative to the modes which enter in the matter dominated era. This fact can 
be used to crudely model the transfer function using the simple form T{k) = 1 for A; < keq and 
T{k) = (keq/k)^ for k > keq, and it depends on just one quantity — the mode which enters the 
horizon at the epoch of matter radiation equality 



In the standard CDM model ^7^0 (= 4.18 x 10^^/i~^) has contributions from the CMBR 
photons and 3 massless neutrino species, and the matter, which is largely made up of CDM 
particles has flmo = ^^cdmo = 1, and we have keq = 0.2/t^Mpc"~^. It should be noted that here, 
and in the rest of the paper, we have ignored the dynamical effect of the baryon density J7bo- 

The decaying neutrino model has all the ingredients of the standard CDM model the only 
difference being that one of the three neutrinos is massive. For neutrino masses rrii, < 1 MeV the 
massive neutrinos decouple when they are relativistic, and after this the comoving number density 
of the neutrinos remains fixed, and at present we expect to find around 112.5 ncutrinos/cm~^ 
(for the present temperature of CMBR Tq = 2.73 K). Once the temperature of the neutrinos 
falls below nii^ the massive neutrinos become nonrelativistic and they contribute to the 
matter density, and, if the neutrinos remain stable up to the present epoch, this contribution 
at present will be ^muo = m^/(93.6h^ eV). So at any epoch before the massive neutrinos 
decay we have oJrijl) = 3.62 x 10~^/i~^ (from the photons and two massless neutrinos) and 
^m{fl) = ot(??)(^^CDMO + ^^mi/o) = + ^mvo), and the epoch of matter radiation equality is at 

teql with 



r?eqi = 2(V2-l) / '\ . (5) 

The universe remains matter dominated from the epoch teqi to t^i when the massive neutrinos 
decay into 2 massless particles referred to as the relativistic decay product. When the massive 
neutrino decays, the matter density falls from cOmii]) = o^(^)(l + ^muo) to ujm{rj) = a{r]), and the 
rest energy of the massive neutrinos gets converted into radiation. This causes ujr{ri) to increase 
from ^IrO to QrO + fld^mi/O (where = a{td)) and the universe becomes radiation dominated for a 
second time. The second radiation dominated era starts at 



1 
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(T+Tl 



1/3 

(6) 



when 

^' + ^r°^'^' [m.Hof (7) 
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and after this the evolution is just hke the standard CDM model except that in this case the 
radiation density is higher. The universe has a second epoch of matter-radiation equality at 



Veq2 



(8) 



and after this the universe is dominated by the CDM particles. The shape of the transfer function 
is determined by the modes which are entering the horizon at the epochs when the universe 
changes from radiation to matter dominated and vice- versa. Expressing these in terms of the 
corresponding mode in the standard CDM model (i.e. keq = 0.2h^ Mpc~^) , we have 



keql = 1.07(1 -I- ^jnuo)keq 



'^eql 
kd 



3.97x10^ sec 

7 (1 + i WOJ 



1/3 



and 



«eq2 



fegq X 



0.866 + 



td 



5.59 X lOiOsec 



2/3 



(1 + QrnuOy^^^^muoh^^^ 



-1/2 



(9) 
(10) 
(11) 

(12) 



For the modes which enter the horizon after the neutrinos decay the transfer function is 
of the same form as the standard CDM transfer function with T{k) = 1 for A; < A;eq2, and 
T{k) = {kcq2/k)'^ for (A;cq2 <k<kd)- Also, as A;eq2 < k^q, the power is reduced on larger scales as 
compared to the CDM model. 

All the modes which enter the horizon in the first matter dominated era grow by the same 
factor and they experience a period of stagnation when the universe is radiation dominated for 
the second time. The amplitude of these modes is suppressed by the factor (r/d/'^cq2)^ relative to 
the modes which enter in the final matter dominated era, and the transfer function has the form 
T{k) = {keq2/kd)^ in the range {kd < k < keqi). The modes which enter the horizon in the first 
radiation dominated era experience a further period of stagnation corresponding to the interval 
between their entering horizon and r^eqi when the universe becomes matter dominated for the first 
time, and for k > keqi the transfer function has the form T{k) = (/ceq2/^d)^(^eqi/^)^- 

In figure 1 we show the contribution to uj(r]) from the various species for a universe with 
h = 0.5, and a decaying neutrino with = 200 eV and td = 10^^ s. In this figure the contribution 
to uj{ri) from the different constituents is plotted against the mode k = Tr/crj which is entering the 
horizon at that epoch. This illustrates how the dominant component changes and which are the 
modes inside the horizon as these changes occur. 

In figure 2 we show the transfer function for this decaying neutrino model based on the 
simple considerations discussed above. In the same figure we have also shown the results of a more 
detailed numerical computation of the transfer function for the same decaying neutrino model. 
In the latter calculation we have numerically followed the evolution of both the metric, and the 
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relevant properties of the different constituents. For the background universe we have solved for 
the scale factor along with the density of the different constituents. For the massive neutrinos 
we have used the distribution function to calculate the background density and this accurately 
follows the evolution of these particles from the relativistic to the nonrelativistic regime. This also 
takes into account the fact that relativistic time dilation causes the faster moving neutrinos to live 
longer in the frame of the cosmological observer. 

We have calculated the evolution of the metric perturbations in the synchronous gauge and we 
have dealt with only the scalar part of the perturbations. For the CDM perturbations we have used 
the pressureless hydrodynamic equations. We have treated the photon-baryon fluid in the tightly 
coupled limit and we have ignored the baryon density in studying the dynamics. For the massless 
neutrinos, the massive neutrinos and the relativistic decay product we have used the collisionless 
Boltzmann equation to follow the evolution of perturbations of the distribution function, and we 
have used this to calculate the perturbation in the density, pressure and the anisotropic stresses. 
This takes into account the fact that the relativistic neutrinos have a very large mean free path 
(of the size of the horizon) and all sub-horizon perturbations in this component are wiped out 
due to the free-streaming of these neutrinos. Sub-horizon perturbations in the massive neutrino 
component grow only after the neutrino becomes nonrelativistic. The evolution of perturbations 
in the other components affects the CDM perturbations only through the metric perturbation and 
the Appendix contains a more detailed description of the equations we have used to follow the 
evolution of the perturbations. 

In figure 2 we also show the numerically computed transfer function for the standard CDM 
model with h = 0.5. We see that the main effect of the decaying neutrinos is to reduce power on 
all scales except on a small range of scales near the first matter radiation equality A;eqi . 

We have normalized the power spectrum for the decaying neutrino models by calculating the 
r.m.s. quadrupole moment Qrms of the angular distribution of the temperature fluctuation in the 
CMBR, and we have normalized the power spectrum to Qrms = IT/iK. This is consistent with 
data from four years of COBE-DMR observations (Bunn & White 1996). On large angular scales 
{6 > 1°) the major contribution to the anisotropy in the CMBR is due to the linear Sachs- Wolfe 
effect (Sachs & Wolfe, 1967) which is an integral of the rj derivative of the metric perturbations 
along the photon's trajectory from the last scattering surface to us. In the standard CDM model 
the photons decouple from the baryons in the matter dominated era, and this integral is restricted 
to the matter dominated era and it can be reduced to a surface term which relates the anisotropy 
in the CMBR to the fluctuations in the gravitational potential at the last scattering surface. In 
the decaying neutrino models we have two radiation dominated era, and there may be signiflcant 
contributions to the CMBR anisotropies from the metric perturbations in the second radiation 
dominated era. For these models it is not possible to analytically calculate the CMBR anisotropy, 
and we have numerically evaluated the integral over the metric perturbations to calculate the 
CMBR anisotropy. We have used this to normalize the power spectrum for all the decaying 
neutrino models and the details are presented in the Appendix. 
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The normalization is significantly different from the standard CDM model in only those cases 
where the second matter-radiation equality occurs quite late i.e. (7?eq2 ~ Vo, where r/o is the value 
of the conformal time at present). The extreme case of such a situation is where the universe 
remains radiation dominated until the present epoch and such cases can be ruled out from age 
constraints. If the universe remained radiation dominated until the present epoch, the age of the 
universe would be 1/2Hq^ as compared to 2/3Hq^ for the matter dominated case, and an age as 
short as ~ lO^'^ Gyr is ruled out by the observations of the oldest globular clusters. We have only 
considered those models where the universe is matter dominated at present. As we shall see in the 
next section, there are models which satisfy this criterion but still have a significant contribution 
from the 'Integrated Sachs- Wolfe effect,' and the normalization of the power spectrum for these 
models is quite different from the standard CDM normalization. 

It is possible to think of the HCDM models as a limiting case of a decaying neutrino model 
with a small neutrino mass and ta ^ to, however, there is a qualitative difference between the 
HCDM models and the decaying neutrino models. For the decaying neutrino models the constraint 
that the age of the universe is required to be c± 2/3Hq^ implies that the energy in the form of the 

relativistic decay products must be negligible at present, and the universe is largely dominated by 
the CDM particles (r^cDM ~ 1) which clump at small scales. This should be contrasted with the 
the HCDM models where 20-30% of the matter at present is in the form of low mass neutrinos 
which do not clump at small scales. 



3. Results 

In this section wc analyze the numerically computed power spectrum for a large class of 
decaying neutrino models and we compare our results with observations at various length scales. 

Allowed models are selected on the basis of comparison to: 

1. the observed cluster abundance (which probes scales ~ 8h~^ Mpc) (White et. al 1993; Eke 
et. al 1996; Henry & Arnaud 1991; Viana & Liddle 1996; Bond & Myers 1996; Pen 1996; 
Borgani et al. 1997; Carlberg et al. 1997), 

2. the peculiar velocity measurements on scales up to 60h~^ Mpc (Bertschinger et al. 1990; 
Courteau et al. 1993; for comprehensive reviews see Strauss &; Willick 1995; Dekel 1994 and 
references therein ), 

3. the three-dimensional power spectrum derived from the APM survey (Baugh & Efstathiou 
1994), which determines the shape of the power spectrum on scales from 1^~^ Mpc to 
50^""^ Mpc, and 
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4. the power spectrum obtained from the Las Campanas Redshift Sm'vey (LCRS) (Lin et al. 
1996). 

We also compare the predictions of the allowed decaying neutrino models with the observed 
abundance of damped Lyman-a clouds at high redshifts. 



3.1. Analysis of the linear power spectra. 



For the scales which enter the horizon after the massive neutrinos decay (i.e. k < kd), the 
massive neutrino model is akin to the standard CDM model, the only difference being the delay 
in the matter-radiation equality due to the additional radiation density contributed by the decay 
products. The delay in matter-radiation equality can be cast as a change in the 'shape parameter' 
r = 5 keqh~^ Mpc for the power spectrum in CDM-like models, and for the decaying neutrino 
models we can similarly use F = 5 kcq2h^^ Mpc, where keq2 can be calculated using equation ( |T2[ ) 
This definition of F is the same as the one used in equation (||) when the value of k^q is expressed 
in terms of QrO and Omo- A point to note is that in the decaying neutrino models F depends on 
the parameters mj, and td through the combination i^(l -|- rriu/ (93.6/i^ eV))mj^ (eq. ^). and there 
is a degeneracy in the relation between m,^, and F. In the limit nii, S> 93.6/i^ eV this relation is 
simpler and we find F depends on just the combination m'^td, which is in agreement with White, 
Gelmini, & Silk (1995), and contrary to the m^td scaling obtained by Bond & Efstathiou (1991). 

A comparison of various numerically computed power spectra with the fit using 
F = 5 kccfih""^ Mpc shows that for large masses this slightly underestimates the value of F, 
and we obtain a better fit using a minor variant of the form for F proposed by White, Gelmini, & 
Silk (1995). 

We find that for large masses rriy ^ 1 keV and small lifetimes td < lOOyr, the power spectrum 
in the range k < l/i~^Mpc~^ is well described by a CDM-like power spectrum 

P{k) = Akx \l + {ak + {hkf/^ + [ckfy] , (13) 
with a = (6.4/F)/i-i Mpc, b = {^/T)h-^ Mpc, and c = {1.7 /T)h'^ Mpc. where our best fit F is: 



F ~ /i X 1 0.15 



, a -|2/3\ -1/2 

rriy Y td 



IkeV/ lyr 



(14) 



which differs from the fit given by White, Gelmini, and Silk (1995) only in a small change of 
the numerical coefficient. The value of normalization A is determined from COBE observations 
and may contain a substantial contribution from the 'integrated Sachs- Wolfe' effect (Section 
(2) and Appendix C). For these decaying neutrino models (rrijy ^ IkeV and td < lOOyr) a 
pure CDM-like fitting function correctly describes the linear power spectrum for the purposes 
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of comparing the theoretical predictions with observations Hke pecuhar velocities of galaxies, 
galaxy-galaxy correlation functions and cluster abundances which probe the power spectrum 
at scales k < l/iMpc~^. At smaller scales {k > l/iMpc~^) these models start exhibiting an 
enhancement of power due to the matter dominated era caused by the massive neutrinos before 
they decay, and a CDM-like fit is inappropriate. 

For smaller neutrino masses rrii, < 1 keV and long lifetimes td ^ lOOyr, the mode kd enters the 
range k < l/iMpc~^a and a CDM-like fit is not adequate to describe all the features of the power 
spectrum at these length-scales. We find that a CDM-like form correctly describes the initial 
deviation (at k ~ kcq2) of the power spectrum from the primordial Harrison Zel'dovich form, but 
it fails on scales ~ kd where there is a significant enhancement of power in the decaying neutrino 
models. The enhancement of power is most pronounced at the scales kd ^ k ^ k^qi which entered 
the horizon when the universe was dominated by the massive neutrinos. The scales k > keqi 
entered the horizon during the era when the massive neutrino was relativistic and the universe was 
experiencing the first radiation dominated phase. The growth of perturbations on these scales is 
suppressed due to free-streaming of the neutrinos and the oscillations in the photon-baryon fluid. 

In figure |3| we show several power spectra for different decaying neutrino models with h = 0.5, 
and my and td such that for all of them the product m^td has the constant value which, according 
to equation (|lj), corresponds to F ~ 0.24. In the range k < 5/iMpc~^ the curve for rriy = lOkeV is 
indistinguishable from a CDM-like spectrum with F = .24. As rriy is reduced the power spectrum 
starts deviating from the pure CDM-like form and it starts exhibiting the various small scale 
features mentioned above and discussed in detail in section 2. The change in small scale power in 
the linear power spectrum can be quantified using an which is the theoretically predicted r.m.s. 
mass fluctuation in randomly placed spheres of radius i?/i~^Mpc. It should be noted that the an 
discussed here is calculated using linear theory and for scales smaller than R = 8 non-linear effects 
have to been taken into account before this can be compared with observations. The change in 
small scale power as the neutrino mass is decreased is shown for various length-scales in figure |^, 
the largest scale being 8/i~^Mpc where the mass fluctuation is observationally very well determined 
and the smallest being .25/i^^Mpc which is comparable to the scales important for the formation 
of damped Lyman-a clouds discussed later in this section. We flnd that as the mass is decreased 
there is nearly a two-fold increase in the small scale power compared to what one expect from a 
pure F = .24 CDM spectrum. At the scale 8/i~^Mpc the mass fluctuation keeps on increasing 
monotonically as rrii, is reduced and this results in some of the low mass models predicting too 
high a value of ag to be compatible with observations. At smaller scales ^ 2/i~^Mpc the mass 
fluctuation initially increases as rrii, is reduced but then starts decreasing again. This happens 
because for low neutrino masses these scales enter the horizon before the massive neutrinos become 
non-relativistic and the growth of fluctuations at these scales is suppressed. 

Figure |^ also shows the APM power spectrum with a bias parameter b = 1.2. Also shown is 
the best flt to LCRS power spectrum (Eq. 23 in Lin et al.'s paper) with b = 1.2. While there are 
uncertainties in the value of the bias parameter for the galaxy surveys, it is believed that the value 
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of as as determined from cluster abundances is relatively free of bias. We have chosen the value 
b = 1.2 for the APM survey to make the observed power spectrum compatible with as = 0.55, and 
we have used this to visually compare the shape of observed power spectrum with the theoretically 
predicted power spectra for various decaying neutrino models. For comparison we have also 
shown the T = 0.5 CDM power spectrum which is obtained in the standard CDM model with 
h = 0.5. While we have shown the linear power spectrum for the various theoretical models, the 
observationally determined power spectrum has significant non-linear effects at A; > 0.2/iMpc~^ 
and a detailed comparison is not possible at these scales until these effects have been taken into 
account. 

We next briefly discuss the COBE normalization of the power spectrum for the different 
models. All the models shown in figure ^ have the same normalization which matches with 
the CDM normalization, and this happens because for all these models the second radiation 
dominated era ends much before the present epoch i.e. (r7eq2 ~ -Ol^). In figure |5| we show the 
power spectrum for the decaying neutrino model with = 100 eV, tj = 5 x lO^^sec for which the 
second radiation dominated phase ends much later (??cq2 ~ -l^o) and this corresponds to F ~ .02. 
We find that for this model the normalization is significantly lower than the CDM normalization 
due to the contribution from the 'integrated Sachs- Wolfe effect'. We have considered a large 
number of different models where T/eq2 ~ -1%) all of which satisfy the constraint that to ~ (2/3)//o 
and we find that in all of them the effect is similar to the case shown in figure ^. 



3.2. Allowed Models 



A very sensitive observation for fixing the value of as is the abundance of rich clusters 
at the present epoch (Evrard 1989). Taking into account various uncertainties due to cluster 
temperatures, cluster x-ray luminosities, and the N-body methods in comparing theoretical 
estimates with observations fix the value of as between 0.5 and 0.9 (White et. al 1993; Eke et. al 
1996; Henry & Arnaud 1991; Viana & Liddle 1996; Bond & Myers 1996; Pen 1996; Borgani et 
al. 1997; Carlberg et al. 1997). Larger values also seem to be ruled out by constraints from the 
pairwise random motions of galactic size dark matter halos (Gelb and Bertschinger 1994). 

For a decaying neutrino model to be consistent with the galaxy surveys (APM, Baugh &; 
Efstathiou 1994; LCRS, Lin et al. 1996) at large scales we have used the broad criterion that 
equation (|l^ ) should predict a value of F in the range .2 < F < .3 for the model, and we have 
discussed a set of models for which F = .24 in some detail. 

In Table 1 we list the value of as for a wide range of decaying neutrino models. A class of 
models allowed by this observation is rriiy = IkeV and td = 100 yr, and all the models obtained 
from this model keeping m^t^ constant for m,^ > 50 eV (figure Q), These models correspond to 
F = 0.24 at large scales. Not only do these models produce a value of as in the right range. 
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but as seen in figure ^, they also reproduce the shape of the observed galaxy power spectrum 
quite well. A comparison of the observed bulk velocities at scales ~ 40 Mpc and ~ 60 Mpc with 
the predictions of these models show that these models may not have enough large scale power 
(figure ^). However, due to systematic uncertainties in the analysis of the data and statistical 
uncertainties due to cosmic variance, the actual peculiar velocities can be much lower than the 
values shown in (figure ^ (Dekel 1994), and therefore a marginal inconsistency of our result with 
peculiar velocity measurements cannot rule out these models. 

For the above class of models (i.e. (mjy/lkeV)^(t(f/lyr) = 100) the value of a% keeps on 
increasing as the mass is lowered below ~ 200eV and and the mass range rriy < 50 eV predicts 
too large a value of as to be compatible with observations of cluster abundances. Models with 
small masses {nii, < 50 eV) are ruled out by the cluster abundance constraints irrespective of the 
value of td- We have run models for < 50 eV with varying over seven decades from 10^ sec 
to 10^^ sec and and we do not find any for which erg < 0.9 The upper limit on the value of td is 
chosen from the consideration of keeping the age of the universe ~ 2/3H^^. The reason why it is 
not possible to get acceptable models in the low mass range is easy to understand. For the decay 
products to substantially delay the final matter radiation equality the massive neutrinos have to 
decay after they have become non-relativistic and they have a significant amount of rest energy to 
pump into the decay products. For this to happen these neutrinos should decay only after they 
have caused the universe to become matter dominated and remain matter dominated for some 
time. The epoch when the first matter radiation equality (?/eqi) occurs depends only on the mass 
of the neutrino (equation ^) and for masses m,y < 50eV the mode fceqi is very close to the mode 
0.2/iMpc~^ which is the scale where the power spectrum is probed by fig. As discussed in section 
2, the power in the range of modes keqi > k > kd is enhanced as a consequence of the first matter 
dominated era, and this gives rise to a large value of cjg. Since feeqi has no dependence on td, the 
value of (Tg too is nearly independent of td for these models. 

We find that by varying both rrii, and td it is possible to construct a large number of models 
all of which predict values of cg which are in the correct range, and some of these are shown in 
Table 1 and figure ||. As one notices, models which give power spectra of very different shapes at 
larger scales can still produce acceptable values of fig. However, most of these models have too 
small a power at large scales, and they do not predict a reasonable value for F. In addition, the 
peculiar velocities predicted by these models are too low to be consistent with observations (figure 
I). 

Going back to the allowed models, we find that an interesting region in the parameter space 
lies in the mass range < 2 keV which, apart from being in consonance with all the observations, 
can give extra power at small scales (figure ^; see also, McNally & Peacock 1996). This will 
result in an early epoch of galaxy formation which might ionize the hydrogen in the intergalactic 
medium, which is seen to be highly ionized up to z ~ 5 (Giallongo et al. 1994). Recent observation 
of high redshift galaxies in the Hubble Deep Field (HDF) indicate that the number density of 
galaxies at z ~ 3 may be comparable to the 2; ~ population (Steidel et al. 1996; Madau et al. 
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1996). Such observations require a high redshift of galaxy formation which occurs naturally within 
the framework of these decaying neutrino models. 

Another possibility with the decaying neutrino scenario is to consider radiatively decaying 
neutrinos. In these models, a small fraction B of massive neutrinos decay into a photon and a 
massless neutrino. It is then possible to directly ionize the intergalactic medium at high redshifts 
by the decay photons (rather than by formation of high redshift objects). There exist a number 
of cosmological and astrophysical constraints on models of radiatively decaying neutrinos. These 
include constraints from the spectrum of the CMBR, the supernova 1987A, the cooling of red 
giants, and the diffuse extra-galactic background of photons, and for a detailed discussion of 
constraints on radiatively decaying neutrinos the reader is referred to Kolb & Turner (1990) and 
references therein. The parameter space for radiatively decaying decaying neutrinos which is 
allowed by all these observations and which is consistent with the observed ionization state of the 
IGM had been studied by Sethi (1997). where it was shown that the acceptable values of B lie 
between 10~^ and 10"''. For this range of B, the allowed values of 1^ and m,y lie between: 

id ~ 3 X lO^^sec (m^/100 eV)-^-^ (15) 

and 

~ 3 X lO^^sec (m^/100eV)"^■^ (16) 

We have calculated the power spectrum for several of the models allowed by the intergalactic 
medium ionization and in figure |^ we show some of the power spectra. The corresponding values 
of (Tg and bulk velocities are shown in Table 1. It is clearly seen that all these model seem to be 
at variance with the observed galaxy power spectrum and peculiar velocity measurements. Figure 
P shows the approximate region of the rrii, — parameter space allowed by the IGM constraints 
(Eq. (|l5|) and (|l^)) as compared to the region of parameter space allowed by observations of 
the large scale structure in the universe, and it is clearly seen that there is no overlap between 
the two regions. For a given mass rriu the lifetimes acceptable for large scale structure formation 
are too short as far as ionizing the IGM is concerned, and as seen in figure these radiatively 
decaying models are ruled out by the shape of the observed power spectrum. However, it is worth 
mentioning that as noted in Sethi (1997), several of these models give acceptable value for 
although they predict an unacceptable shape for the power spectrum. It should be mentioned 
that the gap between the ranges allowed by the structure formation and the IGM ionization is 
too large to be bridged by a change in cosmological parameters. However, the behaviour of the 
two allowed regions as cosmological parameters are varied can be qualitatively understood: for 
instance if the value of h is increased this would, for a given $1^, increase Tgp because Tgp oc 
(see e.g., Miralda-Escude & Ostriker 1990), which means more ionizing photons would be required 
to satisfy the GP tests. If one keeps the value of B, the branching ratio, fixed, even larger values 
of {rriy^td} are needed to ionize the IGM (Sethi 1997). An increase in the value of h would also 
mean more power at small scales (an increase in the value of F), which will have to compensated 
by an increase in {rriu^td}. Therefore, the net results of an increase in the value of h would be to 
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shift up both the allowed regions in figure |8[ Such an effect is difficult to quantify over an entire 
range of {m^jid} because of complicated dependence of the ionizing flux on nii, and t^. 



3.3. High Redshift Objects. 



The abundance of high redshift objects like damped Lyman-a clouds is an important diagnostic 
tool for studying structure formation models (Mo & Miralda Escude 1994; Padmanabhan &: 
Subramanian 1994). The neutral hydrogen column density of damped Lyman-a clouds is 
comparable to the present day spiral galaxies. Also the inferred sizes of these objects suggest that 
the damped Lyman-a clouds are progenitors of the present day spiral galaxies (Wolfe at al. 1992). 
Observations of damped Lyman-a clouds show that a large fraction of baryons in the form of 
neutral hydrogen may already have collapsed into forming these systems at z ~ 3. According to 
Lanzetta, Wolfe, and Turnshek (1995) the quantity ^Im, which is the density of neutral hydrogen 
in damped Lyman-a clouds expressed as a fraction of the critical density, can be fit by a simple 
relation in the redshift range z ~ to z ~ 3 (for go = 0.5): 



Ohi(^) 



0.19 ± 0.04 X W'^h'^ exp(0.83 ± 0.15 x z) . 



(17) 



At higher redshift a decrease in Ohi has been reported (see for instance Storrie-Lombardi et al. 

1995) which suggests that the formation of galaxies commenced around z = 3. For comparison 
with structure formation models one needs to know the mass of damped Lyman-a systems and 
this is highly uncertain. We follow Mo and Miralda-Escude (1994) (see also McNally & Peacock 

1996) in assuming that the minimum mass corresponding to these systems is 10^'^h~^ Mq which 
corresponds to a virial velocity of ~ 50km see" ^. This limit comes from the smallest halos which 
could cool sufficiently rapidly to collapse by z ~ 3 (Efstathiou 1992). Knowing the minimum mass 
one can use the Press-Schechter formalism (Press &: Schechter 1974) to calculate the fraction of 
the total matter which could collapse in forming structures with masses > l(f'^h~^MQ at any 
redshift z. This can be multiplied with fi^ to estimate the density of baryons that has collapsed 
into these objects and this gives 



S^col(^min, Z) =^3 X 1 - erf 



<5c(l + z) 



y(2)a(Mrni„,0) 



(18) 



The only input from structure formation models in equation ( p^ ) is the value of cj(Mjnin,0) 
which is the r.m.s. mass fluctuation in spheres of radius R (~ .3Mpc) corresponding to Mmin 
evaluated at the present epoch. The uppermost curve in figure |^ shows the value of aji at a 
comparable length-scale ( R = .25/iMpc) as a function of the neutrino mass for a class of allowed 
models, and this gives an idea of how this quantity changes for different decaying neutrino models. 
We also use 6c — 1.686 which is the value corresponding to the spherical collapse. It should be 
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pointed out that equation ( [T^ ) gives the collapsed fraction of baryons in all the structures with 
masses > Mmin, though the masses of damped Lyman-a clouds lie between some Mmin and M^ax, 
which probably correspond to the virial velocity ~ 200kmsec~^, comparable to the present day 
galaxies. The collapsed fraction is extremely insensitive to the upper cut off in mass and changes 
negligibly as Mmax is changed from the mass corresponding to a virial velocity of ~ 200kmsec~^ 
to infinity which is the value that equation (|T8| ) assumes. This is expected because the probability 
of forming very high mass objects at high redshifts is exponentially suppressed. 



An important point to note is that while il.Hi in equation (17) refers to the neutral hydrogen 
observed in damped Lyman-a clouds, the quantity being calculated here i.e. Qcoi in equation ( |l8|) 
refers to the total amount of baryons that has collapsed into the damped Lyman-a clouds. It is 
expected that some of the collapsed baryons will be in the form of stars and ionized gas, and 
therefore any acceptable model of structure formation should predict a value for J7coi which is at 
least equal to ilm{z) if not larger. 

In figure ^ we show the collapsed fraction Qcoi for some decaying neutrino models. While the 
baryon density was ignored when calculating the matter power spectrum for these models, we have 
used i^B = 0.05 in equation (|l^). The shaded region in figure ^ shows the observed J7hi (Eq. (p^)) 
for h = 0.5 and the data point at z ~ 4 is taken from Storrie-Lombardi et al. 1995. For comparison 
we also show the collapsed fraction for a HCDM model with ^l^, = 0.25, where is the fraction 
of the matter in the hot component. The power spectrum for the HCDM model was computed 
using COSMICS (Bertschinger &: Bode 1995) with the value i^B = 0.05. As is clearly seen the 
HCDM model seems to be at variance with observations as noted by several earlier authors (Mo & 
Miralda-Escude 1994; Padmanabhan & Subramanian 1994; Ma Sz Bertschinger 1994; McNally & 
Peacock 1996; Ma et al. 1997; Gardner et al. 1997). The two decaying neutrino models considered 
here satisfy all the constraints discussed earlier in this section and the power spectra for these 
models is shown in figure ^. The collapsed fraction for both the models exceeds the observed 
Ohi(-z) which indicates that both these models are compatible with the damped Lyman-a cloud 
observations. Similar conclusions are expected to hold for the other decaying neutrino models that 
pass the various tests discussed earlier this section. It should be pointed out that our computation 
of the matter power spectrum doesn't take into account the baryons. The inclusion of baryons 
reduces power at small scales (see e.g. Hu & Sugiyama 1996) which means that the computed 
fraction shown in figure ^ is an overestimate by 15-30% depending on the value of Qb- However, 
it is evident from figure |^ that a decrement of this level in our computed Qco\ will not affect the 
conclusions of this comparison with observations of damped Lyman-a clouds.. 

Similar constraints on structure formation models can also be derived from the recent 
observation of high redshift galaxies (Steidel et al. 1996; Madau et al. 1996). However, there is 
an even greater uncertainty attached to the masses of these high redshift galaxies and comparison 
with structure formation models may not be so straight forward (Mo & Fukugita 1996) as for 
damped Lyman-a clouds. 



- 18 - 



4. Summary and Discussion 

We have studied large scale structure formation in a 17 = 1 universe with decaying massive 
neutrinos. This variant of the CDM model has two extra parameters — the neutrino mass rriu 
and the lifetime of the massive neutrino td, and by varying these it is possible to introduce extra 
features into the power spectrum. We have computed the power spectrum for a large range of the 
parameters rrij, and t^- Our analysis takes into account the free-streaming of the massive neutrinos, 
and this allows us to study a hitherto unexplored region of parameter space corresponding to low 
neutrino masses and large lifetimes. 

Unlike other models of structure formation in the universe - for example the standard 
CDM model, and its variants like HCDM, ACDM, oCDM - the decaying neutrino model 
allows the possibility of introducing extra features in the power spectrum at both large scales 
{k < 0.01 Mpc~^) and small scales {k > O.lMpc"^). The decay of the neutrino acts to reduce 
the power at large scales (as compared to CDM model) by increasing the radiation content of 
the universe to a value which can be much larger than what is obtained from three relativistic 
neutrino species. This reduces the power at large scales by delaying the matter radiation equality. 
In addition, the late radiation dominated era can affect the overall normalization of the power 
spectrum through the 'integrated Sachs- Wolfe effect' and we find that this can lead to an overall 
reduction of power at all scales. The decaying neutrino models also allow us the possibility of 
enhancing the power over a range of modes at small scales, and this is achieved by varying the 
parameters to suitably adjust the first era of matter radiation equality. 

CDM-like models with T between 0.22 and 0.29 seem to be compatible with most observations 
(Peacock & Dodds 1994). To allow for the uncertainties in comparing the computed linear power 
spectra with the observed non-linear ones, i.e., from APM and LCRS surveys, we allow a larger 
range of F: 0.2 < T < 0.3. We find that it is possible to construct a large number of models 
which predict a value of T in this range and the allowed region of parameter space is shown in 
figure p. We have studied in some detail a class of models whose power spectrum is similar to a 
r = .24 CDM power spectrum at large scales near the turnaround from the primordial Harrison 
Zel'dovich form. We find that for the mass range rrii, > 10 keV, for the entire k range that we have 
studied (i.e. k < 5/iMpc~^), the power spectrum is indistinguishable from a CDM power spectrum 
with r ~ 0.24. As the mass is reduced below lOkeV the power spectrum starts getting extra 
power at small scales and the shape of the power spectrum starts to differ considerably from a 
CDM-like power spectrum with a T fit. We find that the models in the mass range niiy > 50eV are 
roughly consistent with the APM and LCRS power spectrum, and with ag inferred from cluster 
abundances, but with extra power at small scales. We find that the predicted r.m.s. peculiar 
velocities in spheres of radius 40/i~^Mpc and 60/i"^Mpc are somewhat below those indicated by 
observations, but the discrepancy is not as severe so as to conclusively rule out these models. For 
masses < 50 eV the values predicted for erg are too large to be compatible with observations of 
cluster abundances. This conclusion is independent of the choice of t^. and we fail to find any value 
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of trf for which models with rui, < 50eV produces an acceptable power spectrum. 

We also find that for very small masses the neutrinos have to decay very late if this process is 
to have a significant effect on the power spectrum. In these models the power spectrum has too 
much power at small scales and too little power at large scales, and these models are ruled out by 
observations (figures ^ |5| and ^, and Table 1). 

We have compared the predictions of some of the allowed decaying neutrino models with the 
observations of the abundances of damped Lyman-a clouds at high redshifts and we find that the 
decaying neutrino models are compatible with these observations. 

We have addressed the question if any of the radiatively decaying neutrino models which 
can reionize the intergalactic medium is consistent with the observed large scale structure of the 
universe. We fail to find any model which can reionize the IGM and also produce an acceptable 
power spectrum in a 17 = 1 universe. It is still to be seen if such models can be constructed 
by considering open universes or by introducing a cosmological constant. In addition to the 
constraints on radiatively decaying neutrinos considered in this paper, models which predict very 
early reionization are also constrained by CMBR observations at small angular scales (Netterfield 
et al. 1997). However, as we find that these models are ruled out by the structure formation 
considerations presented in this paper, we have not considered the small angle CMBR observations 
here. Incidentally, an early epoch of enhanced galaxy formation may be expected in the allowed 
models with extra small scale power, and this may provide a method of reionizing the IGM. Issues 
related to galaxy formation in a decaying neutrino model are beyond the scope of this paper, 
and they require a more detailed analysis involving N-body simulations and baryonic physics. 
Such an analysis will also be able to put further restrictions on the allowed regions of the my-td 
parameter space. In addition, N-body simulations for the decaying neutrino power spectra will 
permit comparison with the observed power spectrum at scales where the non-linear effects are 
important and this can be used to put restrictions on the allowed models. 

Finally, although we see that it is possible to use the observed distribution of matter in the 
universe to restrict the parameter space of decaying neutrino models, it seems implausible that 
these observation will be able to distinguish a decaying neutrino model from other variants of the 
CDM model in the near future. However, the study of CMBR anisotropics at small angular scales 
holds the promise of being able to discriminate between these models and possibly single out 
the correct one. Two satellite projects planned for the next decade — MAP and Planck surveyor 
(COBRAS-SAMBA) — will map the CMBR sky at a few arc-minutes angular resolutions, and it is 
hoped that these observations will unambiguously validate one of the models and reject the rest. 
In light of this it important to understand the small angular scale anisotropics in the CMBR for 
the various allowed decaying neutrino models, and work is in progress in this direction. 
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A. Formalism 

Here we briefly discuss the equations that govern the evolution of the scale factor and linear 
perturbations in a spatially flat universe which is composed of CDM particles, photons tightly 
coupled to baryons, massless neutrinos and massive neutrinos which have the possibility of 
decaying into some massless (relativistic) decay products. The evolution of both the background 
universe and the perturbations is governed by the Einstein equations 

R^^ = 87rG{Tl^-^SiiT). (Al) 

where the Ricci tensor is calculated from the metric, and the energy momentum tensor Tjf has 
contributions from all the different particle species. In this part of the paper we use c = 1. 

In the synchronous gauge the metric components can be written as 

500 = -a'^irj) goi = gij = a^{ri) [5ij + /ijj(x, r?)] (A2) 

where the zeroth component of the coordinates refers to the conformal time 77 and hij is the metric 
perturbation. Using this in the Einstein equation one obtains the equation for the scale factor and 
the metric perturbation. Using prime to denote derivative with respect to 77, the equation for the 
scale factor can be written as 

a [r^) = H^.f^) . (A3) 

This equation is discussed in some detail in section 2. We next consider the equations for the 
evolution of /t*- the metric perturbation. We proceed by first calculating the perturbation in the 
Ricci tensor caused by the metric perturbation. We only use the components Rq and and 
keeping only terms linear in /i*-, we obtain 

and 

Here h is the trace of the metric perturbation, and we have used the notation j = 
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It is possible to decompose the metric perturbation into scalar, vector and tensor components, 
and at the hnear order these evolve independently. Retaining only the scalar parts, we write the 
metric perturbation as 



/i;(x,r/) 



{2tt)3 3 



/i(k,r?)5j+A(k,r/) K 



(A6) 



where which is the Fourier transform of h corresponds to isotropic dilations or contractions, 
and A, which is the traceless part, corresponds to shear. Expressing the Fourier transform of the 
curvature in terms of these quantities we have 



Ai?0(k,r/) = ^(^/." + (A7) 



and 



Ai?°(k,r?) = ^(M' + A') (A8) 



We also write the Fourier transform of the linear order perturbation to the energy momentum 
tensor as 

ArO(k, 7?) = -Ap(k, r]) Ar](k, r]) = AP(k, rj)S] + ^AQ{k, rj)i5) " 3^) • (A9) 

where Ap(k, ry), AP(k, ry)andA(5(k, rj) are the Fourier transform of the perturbations in the density, 
the pressure and the anisotropic stresses respectively. The perturbations to the components 
also have contributions at linear order, but as discussed below, the final equations that we use do 
not have any explicit reference to ATj*. 

Using the above expressions in the Einstein equation, the q component gives us an equation 
for fi 

+ = -3Hy ( ^ + 3—) . (AlO) 



a \ Pco Pco 

To follow the evolution of A we use the ^ component of the Einstein equation which gives us 



Pco k'^ 



,' + X' = =^%Af^. (All) 



Differentiating this once with respect to r] and using the energy-momentum conservation equation 
Tj".^ = 0, we obtain the following equation for A 

A" + ^ f2A' + p) = 3Hy ( ^ + 3^) . (A12) 
a ^ ' \Pco Pco ) 



Finally, we have the equations ([A3|) , ( [AlO ) and ( A12| ) which we use to follow the evolution of 
a,/i and A respectively. The right hand side of these equations has quantities which refer to the 
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total density, pressure and anisotropic stresses. These have contributions from the different species 
of particles and we have to consider them separately. The total effect is obtained by summing 
over the contribution from the different species e.g. u = ^ LOgpeciesi^P = S ^Pspedes-, etc. In 
the following subsections we separately consider the different kinds of species we have taken into 
account. 



A.l. The ideal fluids. 

Cold dark matter particles and photons which are tightly coupled to baryons can be treated 
as ideal fluids. The energy momentum tensor then has the form 

Tl^ = P5i^ + {P + p)U^'U^ (A13) 

where U^{x,ri) is the bulk 4 velocity of the fluid, and there are no anisotropic stresses i.e. (5 = 0. 

The CDM particles can be considered as dust for which there is no pressure i.e. Pcdm = 0. 
We choose a synchronous coordinate system which moves with the dust particles and hence these 
particles have no peculiar velocities i.e. t^^DM ~ O-I'^^^i^^^)- The energy-momentum conservation 
equation Tcdm o;^ for the CDM component then gives the equation 

<^cdm('7) = a(77)ricDMo (A14) 
for the background density, and the equation 

<5cDM(k,r7) + ^/x'(k,r?) =0 (A15) 

with 

ApcDM _ ^CDM'^CDM 
4 ' 

for the density perturbation. 

The photon-baryon fluid has the equation of state = (l/3)/9^. For this fluid the energy 
momentum conservation equations give the equation 

a;^(r?) = o (A16) 

for the background density. Perturbations in this medium produce bulk flows relative to the 
synchronous coordinate system, and in addition to the perturbation in the density, one has to 
also consider the perturbation to the velocity of the fluid A?7^(a;, 77). Only the divergence of this 
quantity {9^ = J couples to the density perturbation, and the energy-momentum conservation 
gives us the equations 

(5;(k,r7) + ^e^%ri) + ^/.'(k,,?) = (A17) 
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and ^ 

ft,'(k,r/)-^5^(k,r/) =0 (A18) 

which we use to follow the evolution of /S.p^/ p^o = lo^b^ja^ and AP^/ = uJ-^S-y / {?>a^). 



A. 2. Neutrinos. 



It is not possible to treat neutrinos as a perfect fluid and a microscopic description is required. 
Every particle is fully described by its position in space-time and its momentum. Instead of using 
the momentum components in the synchronous coordinate system, it is more convenient to use 
the components of the momentum on the tetrad 

e^ = ^(^r-^/^0^ ^' = a\5l + \hl)dx^ (A19) 

introduced by Bond & Szalay (1983). Here the tetrad index h takes values 0, 1, 2, 3, and the tetrad 
is orthogonal (but not normal). The metric has components ghc = g(efe.ec) = a~'^rjbc on the tetrad, 
and a particle's 4-momentum p has components p = g^e^ on the tetrad. The components g° are 
related by q°'q^r]ab = —a^m'^, and any 3 of the 4 components are sufficient to fully describe the 
momentum state of a particle. We use the 3 spatial components and the zeroth component is 
obtained from these using = \/ q"^ + a?im?, where q^ is used to denote Sijq'^q^ . 

Next, the equation of motion (parallel transport) Vpp = for a particle is used to arrive at 
the equation for the evolution of 

Q^e,(g'^) = -[(Ve,e,)(e'')]gV. (A20) 
Keeping only terms linear in the metric perturbation equation ( A20| ) gives 



q'ebiq') = -^q'q%hl^-h,:) (A21) 

for the spatial components of the momentum g*, and it follows that in the absence of any 
perturbations the 3 spatial components of the momentum remain constant, and they evolve only 
due to the perturbation in the metric. 

The state of a large number of neutrinos of a particular species can be described by a 
distribution function /(x, q, ry) which is the number density of these particles in phase space, and 
its evolution equation 

q-ea{f) + q^ea{q')-^f = (A22) 

follows from the local conservation of particles in phase space. We also consider situations where 
the particles of a particular species decay to produce particles of a different species, and for such 
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situations it is necessary to introduce a source term in equation ( [A22 ). This is discussed later in 
this section. 

The distribution function is next decomposed into two parts, an isotropic function of 
q corresponding to the distribution of particles in the unperturbed universe, and another 
corresponding to the perturbation i.e. 

/(x, q, r,) = f{q, r?) + J/(x, q, r,) . (A23) 

Equation ( |A2^ ) then gives 

|-/(g,7?)=0 i.e. f{q,v) = m (A24) 

for the unperturbed distribution function. For neutrinos the unperturbed distribution function is 
the Fermi-Dirac distribution function 

2 

^"^'^^ = hi [eMq/kB^o) + 1] ^^^^^ 

where the factor 2 takes into account the fact that for every neutrino species there will be 
both particles and anti-particles, hp is the Planck constant, A;^ is the Boltzmann constant and 
T,yO = 2.726°JC/1.4 is the present temperature of the relic cosmic neutrinos. 

For the perturbation, equation ( A22| ) gives us 

Using -F(k, q, rj) to denote the Fourier transform of 5/(x, q, rj) and using a to denote the cosine of 
the angle between q and k, equation ( |A26| ) can be written as 



d p ^ iaqk ^ q 



/i' + (1 - 3a2)A' ^/ = 0. (A27) 



Q-q q^ 6 L*" ' ^ ' ^ dq" 

and we use these equations to follow the evolution of the energy momentum tensor 

n^ = ^ I'-^fd'q. (A28) 



The background density can be written in terms of the distribution function as 

P(^) = 4 / (A29) 



and the various perturbed quantities that appear in the equation for /i and A can be written as 



Ap{k, ri) = ^J q'F{k, q, r,)d\ , (A30) 

AP(k, ^) = ^ / ^^F{K q, V)d\ (A31) 

and ^ 

AQ(k, r/) = ^ - 3a2)%F(k, q, rj)d'q . (A32) 
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A. 2.1. Massless neutrinos. 



For massless neutrinos the calculation is greatly simplified because q^{q,f]) = q which does 
not evolve in time. The integral in equation ( A29| ) is a constant and as a result oJvif]) — ^vo- 
Equation ( A27| ) is solved by the method of characteristics and the solution is 



F;,(k,q,r?) = F,,(k, q, ryj)e 



-iak{ri-rii) 



(A33) 



^(k,^) + (l-3a2)A(k,r/) 



where r]i is the instant when the initial conditions are specified and from which we start following 
the evolution of the perturbations. We use this in equations ( A30| ), ( |A31| ) and ( A32| ) where the 
Sq integral can be done analytically. The angular integrals involve the following relation involving 
Legendre polynomials Pi{x) and spherical Bessel functions ji{x) 



d rl 



2 J- 



e'""^ Pi{a)da , 



(A34) 



and for I > 2 we also define 



1 



8/3 + 12/2 _ 2/ - 3 
+ (6/3 + 3/2-9/)i«_2(x 



(6/3 + 15/2 + 3/ - 6)j/+2(x) - (4/3 + 6/2 + 2l)ji{x) 



il rl 



(A35) 



We use these to obtain 



PcQ PcQ 3a'* 



+ 



V r 



At (k,r?)jo(^(^ -??)) + 2A (k,fi)j2{k{r] - fj)) dfj 



(A36) 



and 



AQ^(k,7?) 

PcO 



'3a4 



{[fi(k,rji)j2{k{r] - r]i)) + A(k, ?7i)J^2(fe(?? - m))] 



+ 



P iKfi)j2ikir] -fj)) + X {k,fi)J^2{k{r] - r?)) dfj 



(A37) 



It should be noted that in obtaining equations ( [A36| ) and ( |A37 ) we have used a particular form for 
the initial perturbation F(k, q, r^j) which has been chosen such that it corresponds to the growing 
mode of the perturbation, and this is discussed in more detail later. 

We use equation (A36) and ( A37|) to follow the evolution of perturbations in the massless 
neutrinos. 
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A. 2. 2. Decaying massive neutrinos. 



For decaying massive neutrinos a source term has to be included in equation ( |A22 ) and we 



have 



oq' td 



(A38) 



where td is the hfetime of the massive neutrinos. The evolution of the unperturbed distribution 
function is governed by 



d 



fd{q, v) 



m. 



■fdiq,r]) ■ 



g2 Qrj"--^- td 
We write the solution of this equation in terms of the variable 



i^{q,v) 



f{q,v) 



■dfj . 



(A39) 



(A40) 



which corresponds to the proper time of a neutrino with spatial momentum q, and it goes over 
to the cosmological time t in the limit q <C aniiy. This variable takes into account the fact that 
the decay of the neutrino is governed by the passage of time in its own rest frame, and not in the 
frame of the cosmological observer. The solution for the unperturbed distribution function is 



fd{q, r]) 



-'<P(q,v)/to 



Uq) 



(A41) 



where fy{q) is the Fermi-Dirac distribution given in equation (A25). We use equation ( A41| ) to 
follow the evolution of the background density of the massive neutrinos and we have 



pdin) = / d^qq°{q,v)fd{q,v) ■ 



The equation for the perturbation is 

-^^fd + -^^fd 
or] q^ 



fq^ 9 - 

-""ijTTJd 



o" 



2q '^dq-"- td<i 
which can be simplified by defining 

5/d(x,q,r?) = e-'^/*''5/d(x,q,r/) 

and 



5fd 



^9d{q,v) = §^f''^i) 

and equation ( [A43| ) can then be written as 



td dq 



i^iq^ff) 



Uq). 



(A42) 

(A43) 

(A44) 
(A45) 

(A46) 
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Using F(i(k, q, r/) to denote the Fourier transform of 5/(i(x, q, 77), the solution of equation 
( A46| ) can be written in Fourier space as 



Frf(k,q,r/) = Fd(k,q,r/i)e 
q n d 



-iakr 



+ 



6Jr,, dq 



9d 



{q,fj) M'(k,r/) + (l-3Q2)A'(k,r/) e-^'^^^-^Uf] (A47) 



where the variable 



-dm 



(A48) 



corresponds to the comoving distance a neutrino travels in the time interval {rj — fj), and we use f 
to denote T{q,fj). For massive neutrinos it is not possible to do the q integral analytically. Doing 
the angular integrals analytically we obtain 



Apd(k, T]) 



27r 
3a4 



d_ 
dq 



Uq) [/x(k,7?i)io(A:r) + 2A(k,r/i)j2(fcr)] 



+ 



-g^jddiq, fi) [a^ (k, fi)jo{k{T - f)) + 2A (k, fi)j2{k{T - f))J dr) 



3APrf(k,ry) = ^ | dg^e"^/*^ Kk, r/,)io(fcr) + 2A(k, r?,)i2(A:r)] 



+ 



m 9q 



and 



9diq,v) /i (k,?7)jo(A:(r-f)) + 2A (k,7?)j2(fc(r-f)) dfj 



3AQrf(k,77) = ^|<ig^e-'^/*^|^/,(g)[/x(k,7?0i2(A:r) + A(k,r?,).^2(A:r)] 



+ 



'r?. dq 



9d{q,v) (k,?/)j2(fc(T - f)) + A (k,??)J^2(A;(T - f)) dr? 



We also use equation ( A47| ) to follow the evolution of 

Anrf(k,r/) = j Fd(k,(i,r])d^q 

and we obtain 



2tt 

T 



dqq e 



dq 



Uq) [fi{k,riMkT) + 2X{k, rii)j2{kT)] 



+ 



-g^9diq, fj) (k, fj)jo{k{T - f)) + 2A (k, fj)j2{k{T - f))J dfj 



(A49) 



(A50) 



(A51) 
(A52) 



(ASS) 



for this quantity which is relevant later for the study of perturbations in the decay product. 

We use equations ( [A49D , ( [A50 ) and ( |A51 ) to follow the evolution of perturbation in the 
massive decaying neutrinos. 
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A. 2. 3. Decay product. 



Consider the decay of a massive neutrino which has momentum q and hence is in motion 
relative to the observers who define the synchronous coordinate system. In its own rest frame the 
neutrino decays by emitting 2 massless particles in exactly opposite directions 1 and —1, and the 
direction 1 is isotropically distributed. We use Qi(q, 1) to denote the tetrad components of the 
momentum of the decay particle which was emitted in the direction 1 in the neutrino's rest frame. 
Although the I's have an isotropic distribution, the Qi(q, l)s will not be isotropically distributed, 
and this happens because of the transformation from the the neutrino's rest frame to the tetrad. 
To avoid confusion with the phase space of the massive neutrino, we use coordinates (x, Q, t]) 
on the phase space of the decay product particles, and we add a source term to equation ( A22| ) 



to follow the evolution of the decay products. This extra term takes into account the fact that 
massive neutrinos with all possible momenta q decay in all possible directions 1 (in the respective 
rest frame) to give rise to 2 decay product particles, and the equation for the evolution of the 
decay product is 

Q«e,(/H(x,Q,r?)) + Q«e,(g*)^/H(x, Q, r/) = 

nil, f ^3 dQ 



td 



J d\^6HQ-QiM)U^,ci,v). (A54) 



For the unperturbed distribution function this gives us 

Q-^fRiQ,v) = ^ I d'q^S\Q-Q,{q,l))Uq,v)- (A55) 
which on integrating over d^Q gives 

^ u;n{v) = ^ fdhUq,v) (A56) 



dv Pcotd 

for the background density of the decay product. For the perturbation we have 

^ / d\^5\Ct-CtM,mU^,^,v)- (A57) 

The initial condition for the decay product is different (i.e. SfR{x,Q,r]) = 0) as we assume that 
the initial density of the decay products is zero. Using this, the solution to equation ( [A57 ) can be 
written in Fourier space as 

Fnik,Q,v) = ^ rT^MQ,r?)L'(k,7?) + (l-3a2)A'(k,7?)le--'=(''-'^)dr/ 
6 Jn OQ L J 



6 Jru dQ 

^ /A^'5'(Q-Qi(q,l)) r a\fi)F,{\.,c^,f,)e-^-^^^--^)df, 

Qtd J 47r Jri, 



(A58) 
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where now a is the cosme of the angle between Q and k, and ^(^(k, q, r/) = e '^/*'*-F(i(k, q, We 
use this to calculate the density perturbation for which we obtain 



ApR(k,r?) 



,AP^(k,ry) 2u^R 



PcO 



3a4 



PcoaHd J " ' 47r 



P (k, fi)joik{r] - fj)) + 2A (k, f])j2{k{r] - 77))J dfj 

(A59) 



where ai(q, 1) is the cosine of the angle between Qi(q, 1) and k. Doing the integral keeping the 
q and 1 dependence of ai(q, 1) is rather complicated. The calculation is greatly simplified if we 
assume that the neutrino is at rest when it decays. Under this assumption the angle ai is the 
cosine of the angle between 1 and k and it no longer depends on q. This assumption is reasonably 
good in situations where the neutrino decays much after it has become nonrelativistic. 

Under this assumption it is possible to analytically do the angular integral d^\, and we obtain 



PcQ 



PcO 

m,, 



3a4 



P (k, r])jo{k{r] - fj)) + 2A (k, f])j2{k{r] - f])) dfj 



PcoaHd 



« iv)joik{v - v))^ndiKv)dfj 



(A60) 



Vi 



and 



PcO 



+ 



3a4 
Pcoa'^td Jr,.. 



P (k, 7))j2{k{r] -??))+ A (k, fi)J^2{k{r] - r/))J drj 
a^{v)h{k{r] - fi))5nd{\<i,fi)dfi . 



(A61) 



where 



(5nrf(k,?7) = / Fd(k,q,r?)d3g. 



The evolution of 5nd(\<i,fj) is governed by equation ( |A53 ) which has been obtained in the previous 
subsection. . We use equations ( A60 ) and ( [A61 ) to follow the evolution of perturbations in the 
relativistic decay product. 



B. Initial conditions. 



Here we briefly discuss the initial conditions for adiabatic perturbations. The initial conditions 
are set at an early epoch when the universe is dominated by the relativistic particles i.e. photons 
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and neutrinos, and the massive neutrino behaves hke a relativistic particle and can also be treated 
as a massless species. The scale factor then evolves as 



a(ry) = r]Ho^n^o + K + 1)^^.0 (B62) 

where Ui, is the number of massless neutrino species, and i^i^Q is the contribution that one massless 
neutrino species would make to the present value of Oq- We also introduce the quantity 



^^07 + {nu + l)f^Oi/ 

which is the ratio of the density of the neutrinos to the total density of the universe in the 
radiation dominated era. 

The initial epoch is also chosen such that all the relevant modes are outside the horizon i.e. 
krji <C 1. In this limit the equations for the evolution of perturbations in the photons and the 
neutrinos are quite simple and can be analytically solved. In the limit krji <C 1 the equation for 
the perturbation in the photon-baryon fluid becomes 

<^;(k,r?)) + ^/x'(k,7?)=0. (B64) 

For each neutrino species we define A/9,y(k, ry) = 3APi,(k, r/) = p^{r])6u(k,r]) and 
3A(5;^(k, r/) = pjy(ry)A^(k, r/), and in the limit kr]i <^ equation ( A27 ) gives us 

i(k,r])) + y{k,i^) = (B65) 



and 



A,(k,r?)) + -A(k,77) =0. (B66) 
15 



We see that the evolution of both S-y and 6i, are governed by the same equation and we can 
combine these two by defining 

The photon-baryon fluid has no anisotropic stresses and there is no contribution from the photons 
to AQ. 

We also have equations ( |A10 ) and ( A12| ), which can now be written as 

" 1 ' 6(5 s 

/u + -/X = ^ (B68) 

r] 

and 

X" + -{2X' +fi') = ^{6 + r^A^). (B69) 
1] f]^ 
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We simultaneously solve equations (|B64[ ), ( |B65| ) , ( [B66| ), (|B68 ) and ( [B69 ) to obtain the analytic 



form of the growing mode in the initial epoch, and the solution for the metric perturbation can be 
written as 

rP'C —If) 
Mk,r/) = -L- and A(k, r/) = ^^^^k, r?) (B70) 

where C is a constant which determines the amplitude of the perturbation at the initial epoch. 
Note that we solve for ^ and A only up to an additive constant, and since we only encounter the 
derivatives of /i and A, the additive constant can be ignored for our purposes. 

Using ( p70| ) we can write the solution for the perturbation in the photon-baryon fluid as 

5^(k,r?)=-^ and 0^(k, r/) = . (B71) 

and for each neutrino species we can write the distribution function as 



6 



F,(k,q,r?) = ^ Mk,r/) + (l-3a2)A(k,r/) ^/,. (B72) 



d 



dq" 



The CDM particles do not contribute to the dynamics in the initial epoch, and they move like 
test particles to which our synchronous coordinate system is attached. We use equation ( |A15 ) to 
obtain ^ 

6cdm{Kv) = -^l^iKv) (B73) 
for the perturbation in the CDM component. 

We use these solutions to fix the initial conditions at the epoch r]i . We also assume that there 
is no significant decay of the massive neutrino prior to the epoch r/j and we set the initial density 
of the decay product to zero. 



C. CMBR anisotropies due to the Sachs- Wolfe effect. 



At angular scales greater than a degree the dominant contribution to the anisotropies in the 
CMBR is largely due to the Sachs- Wolfe effect where the fluctuations in the CMBR temperature 
along any line of sight can be related to the derivative of the metric perturbation /i^j(x, 77) 
integrated along the photons trajectory form the last scattering surface to the observer. For an 
observer located at xq, the fluctuation in the CMBR temperature observed in the direction n is 
given by 

AT I no , 

-7^{n) = -- /i,fe(xo - nfj,fj)n^n''dfi. (C74) 

where r/dec refers to the value of the conformal time at the epoch when the baryons and photons 
decouple, and ijq is the present value of the conformal time. 
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The angular dependence of this temperature fluctuation can be expanded in terms of spherical 



harmonics Y'j™'(n 



AT 

— (n)=^aryr(n). 

Lm 



(C75) 



and the ensemble average of the square of the expansion coefficients gives the angular power 
spectrum 

Q =<| ar|2> (C76) 

which, because the ensemble is statistically isotropic, has no m dependence. Writing the metric 
perturbation in terms of the Fourier expansion and using equation (C74) we write the CMBR 
angular power spectrum as 



47r J (27r)' 



where 



Ai{k,r]) = -2tt daPi{a 



-1 



m 1 



^ (k,r/)e 



iak{rio—fi) 



+ (1 - 3q2)A (k,r/)e 



iak{rio—fi) 



(C77) 



df] . (C78) 



which, using equations ( [A34 ) and ( A35| ), can be written as 



27r /-^o 



Vdc 



H (k,fj)ji{k{r]o -??)) + A (k,?7)J^/(A;(??o - ??)) df] . 



(C79) 



where J^i{k{r]Q — fj) is defined in equation ( A35| ). We use this in equation ( |C77| ) to calculate the 
rms quadrupole Qrms = \/5C2/47r which we use to normalize the power spectrum. 



D. The Numerical Scheme and its Accuracy 



To follow the evolution of the background universe we have numerically solved equation 
(A3) for a{r]), along with equation ( A40 ) for il){q,'q) for a set of values of q, and equation ( A56| ) 
for iOii{rj). The q values have been chosen so that they are appropriate for the Gauss-Laguerre 
quadrature scheme, and we have used this method to evaluate the q integral required to evaluate 
ujdiif) (equation( |A42D ) and the right hand side of equation (A56) at each time step. We have 
used 10 points to evaluate the q integrals and we find that there is no significant improvement if 
we increase the number of points. We have used an adaptive step-size fifth order Runge-Kutta 
subroutine 'odeint' (Press et. al. 1992) to follow the time evolution of the set of coupled ordinary 
differential equations. Along with the above mentioned quantities we have also evolved the 
quantities T{q,rf) and ^'^'{l^f]) which are required to follow the perturbations. The intermediate 
values of all these quantities are stored. We fit them with a cubic spline and use the intermediate 
values in studying the evolution of the perturbation. 
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The time-steps have been chosen so as to achieve a relative accuracy of 10-"^. For all the 
models the background universe is first evolved using f^cDMO = 1- It sometimes happens that 
present contribution from the decay products is quite large and the total becomes greater than 
one. In such cases we reduce the value of Ocdmo and evolve the background again, and we keep 
on iterating the process until the value of JIq converges to within 1 it 0.001. 

For the metric perturbation we numerically solve the two second order differential equations 
( AlOj ) and ( A12| ) by converting them into four first order equations. These equations are solved 



together with equation ( |A15| ) for the CDM perturbations, and equations ( A17| ) and ( |A18| ) for the 



perturbations in the photon-baryon fluid. This system of 7 differential equations is evolved using 
'odeint' and the intermediate values of /x and A at each time step are recorded and we fit these 
by a set of overlapping cubic polynomials and these are used in following the perturbations in the 
neutrinos. The metric perturbations are coupled to the perturbations in the neutrinos, and these 
have to be evaluated separately at each time step. For the massless neutrinos we numerically 
evaluate the integrals in equations ( A36| ) and ( |A37 ) for every time step in 'odeint'. Similarly, for 



the massive neutrinos we numerically evaluate the fj integrals in equations (A49), ( A50| ), (A51) 



and ( A53D for a set of values of q and we do the q integrals by a Gauss-Laguerre quadrature. For 



the decay products we have numerically evaluated the integrals in equations (|A60D and ( A61 ) at 
every time-step in 'odeint'. 

We have used 10 points for the q integrals, and we find that increasing the number of points 
does not significantly change the results for the feasible models in the range of k that we have 
considered. 

The results of the computations described above yield the matter transfer function 
T{k) cx| (5(k, ryo) p. We multiply this with k — the primordial Harrison Zel'dovich spectrum to 
obtain the power spectrum. This is normalized using equations ( p77 ) and (|C79 ) which we use to 
calculate Qrvtis • 

In our analysis we have ignored the baryon density and we use = 0. In addition, when 
following the evolution of the dark matter perturbations we have treated the photons as being 
tightly coupled to baryons until the present epoch. Finally, we have altogether ignored the 
interaction of the photons with the electrons and baryons when calculating the CMBR anisotropy. 
All these assumptions cause a few percent error in our results. These effects have been studied in 
detail for several models (see for instance Bond 1996), and it is found that the effect of baryons 
can be included by scaling the value of F as F x exp(— 21^^). Primordial nucleosynthesis imposes 
the restriction 17^^^ — 0.01, and for h = 0.5 one expects a 8% error if the baryons are left out. 

Treating photons as tightly coupled to baryons causes the sub-horizon scales perturbations 
in photons to continue to oscillate with the same amplitude even after recombination. In reality, 
once the photons decouple from the baryons their mean free path becomes comparable to the size 
of the horizon and sub-horizon perturbations in the photons start getting wiped out as a result of 
the free-streaming. However, the photons decouple in the matter dominated era where they play 
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no role in the dynamics of the dark matter perturbations, and these effects are negligible. 

To check the accuracy of our numerical code we have compared our transfer function for 
several models with the runs of COSMICS (Bertschinger and Bode 1995) in the limit ^ 0- 
For CDM models we get an agreement to better than 5% for k ^ 1 Mpc~^. A comparison of the 
transfer functions for various HDM and HCDM models allows us to test the reliability of our 
treatment of the massive neutrinos and we find that for the HDM model, in the range k < 3Mpc~^ 
our transfer function differs from the transfer function calculated using COSMICS by less than 2%, 
which indicates that much of the error in our analysis come from the assumption of tight coupling. 
This results in larger errors in HCDM models, for which, in the range O.lMpc"^ < k < lMpc~^, 
the error is within 10%. 

We have not pinpointed the exact cause of this discrepancy, but some of the possible sources 
are discussed below. One possible source of the error could be the fact that in dealing with 
the massive neutrinos one has to consider the evolution of neutrinos with different momentum 
separately and then integrate over the momentum. We have used only 10 values of momentum 
and the integration was done using a Gauss-Laguerre quadrature where a few of the points 
come with very low weights. COSMICS performs this integration using 8th-order Newton-Cotes 
method with as many as 128 points. This is one of the sources of the error and it may be possible 
to overcome this by using some other quadrature scheme (eg. Bond & Szalay, 1983) and by 
using more points. We have tried doubling the number of points but it does not make a very 
big difference in the results in the range k < 0.5Mpc~^. Another possible cause for difference 
could be the fact that we have done all numerical integrations at a relative accuracy of 10~^ as 
compared to the relative accuracy of 10~^ in COSMICS. A point to be noted is that COSMICS 
uses various moments of the coUisionless Boltzmann equation in order to follow the evolution of 
the perturbation in the neutrinos, and this involves a truncation which is externally enforced. 
Our method does not involve such a truncation as we use the analytic solution of the coUisionless 
Boltzmann equation which is based on the method of characteristics, but it has an extra cost as we 
have to do an integration over the entire past at every time step. Our treatment also has another 
added advantage in that it involves only differential equations, and does not involve any algebraic 
equations, and we would expect it to be more stable compared to methods based on a combination 
of algebraic and differential equations. A little experimentation with the initial conditions shows 
that when they are set arbitrarily (i.e. a mixture of the growing and decaying modes) the solution 
goes over to the growing mode as expected and the effects of the decaying mode die away, showing 
that our numerical scheme is indeed stable and does not dependent critically upon any fine tuning 
of the initial conditions. 

Finally we note that a 10% inaccuracy is acceptable in calculating the matter transfer 
function as it is well below the observational uncertainties. More accurate computations of 
decaying neutrino transfer functions may be required in the future as observations become more 
accurate. Also, a more intensive treatment of the massive neutrinos is required at smaller scales 
which will be essential to address issues related to galaxy formation in the decaying neutrino 



models. A need for more accurate computing will also arise when computing CMBR anisotropies 
for comparing with proposed future observations at small angular scales. Work is currently in 
progress at improving the present numerical scheme so as to be able to address these questions. 
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Fig. 1. — The contribution to uj{r]) from the various components is shown as a function of the mode 
k = tt/ {crj) which enters the horizon at the epoch 77. 
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Fig. 2. — The sohd curves show the transfer function for a decaying neutrino model with h = 0.5, 
777,1/ = 200 eV and ta = 10^^ s. The smooth curve is the result of the numerical computation whereas 
the jagged curve is based on the crude approximation to the transfer function discussed in section 
2. The dashed curve shows the transfer function for the h = 0.5 CDM model. 
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Fig. 3. — This shows the power spectra for several decaying neutrino models all obtained by varying 
m and keeping m^tfi a constant. The CDM power spectrum is shown for comparison. Also shown 
is the APM power spectrum (filled squares) and the best fit to the LCRS power spectrum (dashed- 
dotted thick curve). 6 = 1.2 has been assumed in plotting the APM and LCRS power spectra and 
h = 0.5 is used for all the power spectra shown here. 
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Fig. 4. — This shows how cjr {R in h ^ Mpc) changes as a function of for a class of models for 
which m^(keV)id(yr) = 100. The corresponding values for the standard CDM model are: ag = 1.26, 
(74 = 2.32, (71 = 3.81, C70.5 = 5.75, and c7o.25 = 10-89. 
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Fig. 5. — This shows the power spectrum for several models which give acceptable value of as 
(Table 1). Most of these models are ruled out on comparison to the APM and LCRS power spectra 
which have been plotted here for b = 1.2. These models also fail to make reasonable predictions 
for the peculiar velocities (Table 1). h=.5 has been used here. 
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Fig. 6. — Here we show the r.m.s. bulk velocity in spheres of radius R after smoothing with a 
Gaussian of 12h~^ Mpc. This is plotted as a function of R for the following decaying neutrino 
models, all of which use h = .5: (1) nii, = 1 keV, = 3.15 x 10^ sec {solid curve), (2) nii, = 200 eV, 
td = 10^^ sec {dotted curve), (3) = 150 eV, ta = lO-'^^sec {short dashed curve), and (4) 
777-1/ = 100 eV, td = 5 X 10^^ sec {long dashed curve). The data points shown, V40 = 388 it 67km/sec 
and Veo = 327 ± 88, are taken from Bertschinger et al. (1990) 
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Fig. 7. — Here we show the power spectrum for some of the radiatively decaying neutrino models 
(where a smah fraction of neutrinos B ~ 10^^-10^'' decay into photons) which satisfy various 
observations of the ionization of the intergalactic medium at high redshifts. 
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Fig. 8. — The parameter space aUowed by structure formation bounds (shaded region) corresponding 
to 0.2 < r < 0.3 where r(m,y,trf) is given by Eq. (14) is shown along with the region on the rriy-td 
plane allowed by IGM considerations (i.e. Eqs. (15) and ([T6|)) (hatched region). 
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Fig. 9. — The density of collapsed baryons f^coi for several structure formation models is compared 
to the density of neutral hydrogen Ohi observed in the damped Lyman-a clouds. The shaded region 
corresponds to Eq. (|l^) ( Lanzetta, Wolfe, and Turnshek 1995). The data point at z ~ 4 is taken 
from Storrie-Lombardi et al. (1995). 




Table 1: Values of as and r.m.s. bulk velocity in km/s in spheres of radius 40/i ^ Mpc and 
60/i~i Mpc after smoothing with a Gaussian of 12h~^ Mpc. 
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